function [Z,W,Wvar,invtausq_w,theta,theta_var,invtausq_theta,ELBO] = ...
        run_modQTL_association(EXPR,GENO,ANNOT,CIS,OuterIter,InnerIter)
%  [Z,W,Wvar,InvTauSqW,Theta,ThetaVar,InvTauSqTheta,ELBO] = run_modQTL_association(EXPR,GENO,ANNOT,CIS,OuterIter,InnerIter)        
%
% EXPR      = gene x individual expression
% GENO      = snp x individual genotype
% ANNOT     = snp x K prior SNP annotation
% CIS       = gene x snp cis-relationship
% OuterIter = optimization steps of outer loop
% InnerIter = optimization steps of inner loop
%
% Z         = snp x samples conditioned PIP
% W         = gene x snp association
% Wvar      = variance of W
% THETA     = sample x (K+1) prior log-odds parameters
% ELBO      = sample x 1 (log-importance)
%
% code: Yongjin Park, ypp@csail.mit.edu
%

    % ****************************************************************
    [Ngene,Nind] = size(EXPR);
    [Nsnp,~] = size(GENO);

    ANNOT = [ones(Nsnp,1),ANNOT];
    K = size(ANNOT,2);

    TOL = 1e-2;
    SIGSQ = var(EXPR(:));

    CIS = logical(CIS);

    % ================================================================
    % common starting point
    prior = -1 * ones(Nsnp,1);

    [W,Wvar,Z,invtausq_w,llik,elbo] = modQTL_association_data(EXPR,GENO,CIS,prior,SIGSQ,10,InnerIter);

    fprintf(1,'Iter = 000, LLIK = %.4e, ELBO = %.4e, Pr(Z) = %.2e\n', llik(end), elbo, mean(Z));
    elbo_trace = NaN(OuterIter,1);

    for outer = 1:OuterIter,
        % recompute prior odds ratio
        [theta,theta_var,invtausq_theta,prior,llik] = modQTL_association_prior(Z,ANNOT,10,InnerIter);
        % recompute data model
        [W,Wvar,Z,invtausq_w,llik,elbo] = modQTL_association_data(EXPR,GENO,CIS,prior,SIGSQ,0,InnerIter,W,Z,invtausq_w);
        elbo_trace(outer) = elbo;

        fprintf(1,'Iter = %03d, LLIK = %.4e, ELBO = %.4e, Pr(Z) = %.2e\n', outer, llik(end), elbo, mean(Z));

        if outer > 5 && abs(mean(elbo_trace((outer-5):(outer-1))) - elbo_old)/abs(1e-10-elbo_old) < TOL,
            fprintf(1,'Converged\n');
            break;
        end

        elbo_old = elbo;
    end

    ELBO = elbo_trace;

end
